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Abstract 

We study Lissajous curves in the 3-cube, that generate algebraic cu- 
bature formulas on a special family of rank-1 Chebysliev lattices. These 
formulas are used to construct trivariate hyperinterpolation polynomi¬ 
als via a single 1-d Fast Chebysliev Transform (by the Chebfun pack¬ 
age), and to compute discrete extremal sets of Fekete and Leja type 
for trivariate polynomial interpolation. Applications could arise in the 
framework of Lissajous sampling for MPI (Magnetic Particle Imaging). 
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1 Introduction 

During the last decade, a new family of points for bivariate polynomial in¬ 
terpolation has been proposed and extensively studied, namely the so-called 
“Padua points” of the square; cf. @ii5i Buna- They are the first known op¬ 
timal nodal set for total-degree multivariate polynomial interpolation, with 
a Lebesgue constant increasing like (logn) 2 , n being the polynomial degree. 

One of the key features of the Padua points, that has allowed to construct 
the interpolation formula, is that they lie on a suitable Lissajous curve, such 
that the integral of any polynomial of degree 2 n along the curve is equal 
to the 2d-integral with respect to the product Chebyshev measure. More 
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specifically, the Padua points are side contacts and self-intersections of the 
Lissajous curve. 

Motivated by that construction, in the present paper we try to extend 
the Lissajous curve technique in dimension 3. Since the resulting curve is 
not self-intersecting, we cannot obtain total-degree polynomial interpolation. 
On the other hand, we are able to generate an algebraic cubature formula 
for the product Chebyshev measure, whose nodes lie on the Lissajous curve 
thus forming a rank-1 Chebyshev lattice (on Chebyshev lattices cf., e.g., 

MY 

By such formula we can perform polynomial hyperinterpolation, which is 
a discretized orthogonal polynomial expansion m, and can be constructed 
by a single 1-dimensional Fast Chebyshev Transform along the curve. More¬ 
over, since the underlying Chebyshev lattices turn out to be Weakly Ad¬ 
missible Mehes for total-degree polynomials (cf. f7]), we can extract from 
them suitable discrete extremal sets of Fekete and Leja type for polynomial 
interpolation (cf. [6]). We provide a Matlab implementation of the hyperin¬ 
terpolation and interpolation scheme, and show some numerical examples. 
Applications could arise within the emerging field of MPI (Magnetic Particle 
Imaging), cf. [25] . 


2 3d Lissajous curves and Chebyshev lattices 


Below, we shall denote the product Chebyshev measure in [—1, l] 3 by 

1 


dX = w(x)dx , vj(x) = 


Vi 1 - ~ x‘i)(l - x‘i) 


(1) 


Moreover, P| will denote the space of trivariate polynomials of degree not 
exceeding k, whose dimension is dim(P^) = (k + l)(k + 2 )(k + 3)/6. 

Following the lines of the construction of the Padua points, the strategy 
adopted is to seek a Lissajous curve such that the integral of a polynomial 
in P| n with respect to the Chebyshev measure dX is equal (up to a constant 
factor) to the integral of the polynomial along the curve. To this purpose, 
the following integer arithmetic result plays a key role. 


Theorem 1 Let be n € N + and ( a n ,b n ,c n ) be the integer triple 


[P"ni • C n ) 


(fn 2 + \n, fn 2 + n, |n 2 + |n + l) , n even 
(|n 2 + |, |n 2 + |n - |n 2 + |n + |) , n odd 


( 2 ) 


Then, for ever integer triple ( i,j , k), not all 0, with i,j, k > 0 and i+j + k < 
m n = 2 n, we have the property that ia n / jb n + kc n , jb n / ia n + kc n , 
kc n ^ ia n + jb n . Moreover, m n is maximal, in the sense that there exist a 
triple ( i*,j*, k*), i* + j* + k* = 2n + 1, that does not satisfy the property. 
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Proof. See the Appendix. 


Proposition 1 Consider the Lissajous curves in [—1, l] 3 defined by 

l n (0) = (cos(a n 9),cos(b n 9),cos(b n 9)) , 9 € [0,7r] , 

where ( a n ,b n ,c n ) is the sequence of integer triples (U|). 

Then, for every total-degree polynomial p £ P 3 n 


[ p{x) w{x)dx = 7T 2 [ p{i n {9))d9. 

' [— 1 , 1 ] 3 Jo 


(3) 

(4) 


Proof. It is sufficient to prove the identity for a polynomial basis. Take 
the total-degree product Chebyshev basis Ti{x\)Tj(x 2 )T k {xfi), i,j,k > 0, 
i + j + k < 2n. For i = j = k = 0, 0 is clearly true. For i + j + k > 0, by 
orthogonality of the basis 


'[-id] 3 


T i {xi)T j {x 2 )T k {xfi)w{x)dx = 0 . 


On the other hand, 


f 


Ti(cos{a n 9))Tj(cos(b n 9))T k (cos(c n 9)) d9 


= / cos (ia n 9) cos (jb n 9) cos (kc n 9) d9 
Jo 


sin {{ia n - jb n - kc n )9) 


iO'n jbn kc n 
sin((m n - jb n + kc n )9) 


sin((ia n + jb n - kc n )9) 


i^n jbn T kc n 


+ 


0 + ia n +jb n -kc n 
sin((m r) , + jb n + kc n )9) 


g iO"n H - jb n T kc n 

Now, the fourth summand on the right-hand side is zero since ia n + 
jb n + kc n > 0, and thus the whole right-hand side is zero if (and only if) 
ia n - jb n - kc n 0, ia n + jb n - kc n 0, ia n - jb n + kc n 0, which is true 
by Theorem 1 since i + j + k < 2n. □ 


Corollary 1 Let be p € ^n{9) the Lissajous curve m and 

{ | n 3 + | n 2 + n , n even 

|n 3 + |n 2 + |n , n odd 

Then 

u 

p(x)w(x)dx = 'Y^w s p(£ n (9 s )) , 
s=o 



(5) 

( 6 ) 
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where 


s = 0,... ,£t , 


(7) 


Wa = vr cj s 


with 

H = v , 0s 
or alternatively 


(2s + l)-7r 
2^ + 2 




7r 

/i + 1 


s = 0,, 


// = i/ + 1 


— , s = 0,... ,/z , 
h 


i r 


7T 


Wo = w M = — , UJ S = - , s = 1,..., At - 1 . 


2/z 


M 


( 8 ) 


(9) 


Proof. Observe that by Proposition 1 and the change of variables t = cos(0) 
f p(x) w(x)dx = it 2 f p(£ n (9))dd 

J [ —1, l] 3 Jo 

/ ! rjf 

^p(T an (t),T bn (t),T Cn (t))^==, 

where p(T an (t),Ti }ri (t),T Cn (t)) is a polynomial of degree not exceeding 

2v = ma x{ia n +jb n + kc n , i,j,k > 0 , i+j + k < 2n} = 2 n max{a„, b n , c n } . 

The conclusion follows by using the classical Gauss-Chebyshev or Gauss- 
Chebyshev-Lobatto univariate quadrature rules, cf. ((HJ) and (J9]) respectively, 
which are exact up to degree 2v + l using the n + 1 nodes t s = cos (9 S ) and 
weights oj s , cf., e.g., [261 Ch. 8]. □ 

Remark 1 (Chebyshev lattices). We observe that {t n (9 s )\. s = 0, 
are 3-dimensional rank-1 Chebyshev lattices (for cubature degree of exact¬ 
ness 2 n) in the terminology of TJ]. As opposite to IS 1 , where Chebsyhev 
lattices are generated heuristically by a search algorithm, here we have a 
formula to generate rank-1 Chebyshev lattices for any degree. 

2.1 Optimal Tuples and Homogeneous Diophantine Equa¬ 
tions 

An algebraic trivariate polynomial of degree N restricted to the Lissajou 
curve £ n (9) is a trigonometric polynomial of degree lVmax{a n , b n , c n } = 
Nc n . Hence it is some interest to have a triple for which the conclusion 
of Theorem 1 holds with its maximum as small as possible. Indeed, we 
conjecture that the triples ([21) are optimal in this sense. 
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Conjecture 1 Suppose that (a, b , c) is a triple of strictly positive integers 
such that max{a, b, c} < c n , with c n given by ©■ Then there exists an 
integer triple ( i,j,k ), not all 0, and i + j + k < 2 n, such that either ia = 
jb+kc, jb = ia+kc, or kc = ia+jb. In other words, the triples © are those 
satisfying the conclusion of Theorem 1 having the minimum maximum. 


We do not as yet have a proof of this conjecture but can provide a lower 
bound for the minimum maximum of such “good” triples with the correct 
order of growth in n. 

First observe that the conditions of the conclusion of Theorem 1 may 
be expressed more succinctly in terms of a homogeneous linear Diophantine 
equation. 

Lemma 1 Suppose that (a, b, c) is a triple of strictly positive integers. Then 
there exists an integer triple ( i,j , k), not all 0, and i + j + k < N, such that 
either ia = jb + kc, jb = ia + kc, or kc = ia + jb iff there exists an integer 
triple (x, y, z ) € Z 3 such that |x| + \y\ + \z\ < N and xa + yb + zc = 0. 

Proof. If, for example, ia = jb + kc, then — ia + jb + kc = 0 and we may 
take x = —i, y = j and z = k. On the other hand, if xa + yb + zc = 0 then 
not all of x, y and z can have the same sign. There being an odd number of 
them, two of them have the same sign and the other the opposite sign. By 
multiplying by —1 if necessary, we assume that the single sign is negative. 
For example, if it is x that is negative, we may write — xa = yb + zc and 
take i = —x, j = y and k = z. □ 

The classical Siegel’s Lemma (see e.g. [Ml p. 168]) gives a bound on 
the order of growth of “small” solutions of homogeneous linear diophantine 
equations. We may adapt this to our situation to prove 


Lemma 2 (A version of Siegel’s Lemma) Suppose that 1 < n € Z+. Sup¬ 
pose further that a = [a±,a 2 , ■ ■ ■,a d \ € with m > 0, 1 < i < d, is such 
that 

max{a} < M 


where 


H n * d )\ 


2 (=0(n d -'j). 


Then there exists 0 / x € 7L d such that l x »l — 2 n and 


d 

^2Xi a i = 0 . 

2=1 


Proof. Let Sd C denote the set of positive tuples 0 / z £ Z^ such that 
Ef=i z i < n. Then #(S d ) = ( nJ f) ~ 1- 
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Consider the map F : if —>• Z given by 


d 

F(z) := y^ajZj. 

i= 1 


Then F(Sd ) C [1 ,nM] and hence 

#(F(S d )) < nM. 


But 


nM = n 


1 

fn + d\ 

n 

K n ). 



< 


fn + d\ 

V d ) 


2 n < 


CT) 


i, 


i.e., 


#(F(S d )) < #(5 d ). 


It follows from the Pigeon Hole Principle that there exists two different 
tuples t/ 1 ) f € S d such that 


F(yV)=F(yW), 


i.e., 




i =1 

The tuple x := y^ — y ^ has the desired properties. □ 

In our context it means that the minimum maximum of “good” tuples 
is at least 




2 (=0(n'<- 1 )). 


3 Hyperinterpolation on Lissajous curves 

We shall adopt the following notation. We denote the total-degree orthonor¬ 
mal basis of Pn([— 1) l] 3 ) with respect to the Chebyshev product measure 
© by 

= f i (x 1 )f j (x 2 )f k (x 3 ) , i,j,k> 0, i + j + k < n , (10) 

where T m f) is the normalized Chebyshev polynomial of degree m 


T m (-) = a m cos(m arccos(-)) , a m 
with the convention that sign(0) = 0. 


1 + sign(m) 


7T 


m > 0 


( 11 ) 
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We recall that hyperinterpolation is a discretized expansion of a func¬ 
tion in series of orthogonal polynomials up to total-degree n on a given 
d-dimensional compact region K , where the Fourier-like coefficients are com¬ 
puted by a cubature formula exact on P 2 n (K). It was proposed by Sloan 
in the seminal paper m in order to bypass the intrinsic difficulties of poly¬ 
nomial interpolation in the multivariate setting, and since then has been 
successfully used in several instances, for example on the sphere |24| . 

Given a function / € C([—1,1] 3 ), in view of the algebraic cubature 
formula ©, the hyperinterpolation polynomial of / is 


Unfix) = ^ Ci,j,k^i,j,k( x ) > (12) 

0 <i-\-j-\-k<n 

where 

A 4 

C it j,k = 22 W s f{Zn{0s)) 4>ij,k(£n(0s)) ■ (13) 

s =0 

Observe that by construction T~L n f = / for every / € P;[, i.e., H n is a 
projection operator. Among the properties of the hyperinterpolation oper¬ 
ator, not depending on the specific cubature formula provided it is exact 
up to degree 2n for the product Chebyshev measure, we recall the following 
bound for the L 2 error, 


\\f-Unf\\2<27T 3 E n (f), E n (f)= inf 11/-Plloo . (14) 

pSPn 


Consider the uniform operator norm (i.e., the Lebesgue constant) 


H'Hnll = sup 
i¥0 


WnfWoo 

ll/lloo 


max Y ^w s \K n (x,£ n (9 s ))\ , 
*€[-!,i ] 3 


(15) 


where K n (x,y) = J2o<i+j+k<n^i,jA x )^i,jMy) is tlie reproducing kernel of 
P^ with respect to the product Chebyshev measure (JTJ) , cf. [22]. 

In [T7] the bound \\1-L n \\ = 0((y/n) 3 ) has been obtained, as a consequence 
of a general result connecting multivariate Christoffel functions and hyper¬ 
interpolation operator norms. On the other hand, by proving a conjecture 
stated in m , the fine bound 

||Wn|| =C>((logn) 3 ) (16) 

has been provided in [35]J, which corresponds to the minimal growth of a 
polynomial projection operator, in view of [32]. Since "H n is a projection, 
we get the L°° error bound 

11/-^n/||oo = 0((logn) 3 E n {f)) • (17) 

We show now that the hyperinterpolation coefficients {Cjjfc} can be 
computed by a single 1-dimensional discrete Chebyshev transform along the 
Lissajous curve. 
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Proposition 2 Let be f € C{[— 1,1] 3 ), ( a n ,b n ,c n ) the sequence of integer 
triples HjJ ; and v, n, {0 S }, oj s , {u; s } as in Corollary 1. The hyperinterpola¬ 
tion coefficients of f generated by fW can be computed as 


IT 


_ " _ _ _ I 7«1 | 102 , 103 . /a4 

— . CTia n CTjb n Crkc n I + + + 


7«2 I 7q? 3 , 7c*4 


cr, 


Q1 


O', 


a 2 




<23 






(18) 


O'l — jb n -\- kc n , cx 2 — -|- jb n kc n | , 

O3 — |*®n J ^n| “I - ^C n , O4 — ||^On J^n| kc n | , 

where { 7 m } are tde first v + 1 coefficients of the discretized Chebyshev ex¬ 
pansion of f(T an (t),T bn (t),T Cn (t)), t G [-1,1], namely 

7m = ^w s f m (T s ) /(Ta n (r s ),r fen (r s ),r Cri (r s )) , (19) 

s=0 


m = 0 , 1 ,... , 17 with t s = cos( 0 s ), s = 0 , 1 ,... , /i. 


Proof. By the change of variables 6 = arccos(f) which gives 
Ue) = (T an (t),T bn (t),T Cn (t)), 

and by the classical identity T h (t)T k (t) = ~ ( T h+k (t ) + T\ h _ k \(t)) (cf., e.g., 
[26 . §2.4.3]), we get 

f{&n{9))<i>i,j,k{£n{6)) = /(^(arccos(t))) f ian {t)f jbn (t)f kCn (t) 

= /(^n(arccos(t))) a ian a jbn a kCn i (T ai (t) + T a2 (t) + T as (t) + T a4 (i)) , 

and hence we have (USD-CD in view of da, and the fact that { 77 } are 
the nodes of the univariate Gauss-Chebyshev or Gauss-Chebyshev-Lobatto 
formula, with weights oj s , cf. (f8|)-([9|). □ 

Remark 2 (Lissajous sampling). Hyperinterpolation polynomials on d- 
dimensional cubes can be constructed by other cubature formulas for the 
product Chebyshev measure, that can be more efficient in terms of number 
of function evaluations required at a given exactness degree. For example, 
a formula of exactness degree 2 n with 0(n 4 ) nodes for the 3-cube has been 
provided in [20], and used in a FFT-based implementation of hyperinterpo¬ 
lation. Other formulas, in particular Godzina’s blending formulas [23], that 
have the lowest cardinality known in d-dimensional cubes, have been used in 
the package m- All such formulas are based on Chebyshev lattices of rank 
greater than 1, that are suitable unions of product Chebyshev subgrids. 


A first advantage of rank-1 Chebyshev lattices, as observed in general 
in [T3], is that a single 1-dimensional FFT is needed to compute the hyper¬ 
interpolation polynomials. In the present context of sampling on Lissajous 
curves of the 3-cube, this is manifest in Proposition 2. 

On the other hand, one of the most most interesting features of hper- 
interpolation on Lissajous curves arises in connection with medical imaging 
applications, in particular with the emerging 3d MPI (Magnetic Particle 
Imaging) technology. Indeed, Lissajous sampling is one of the most com¬ 
mon sampling methods within this technology, since it can be generated by 
suitable electromagnetic fields with different frequencies in the components, 
cf., e.g., [25], [28]. Choosing the frequencies (J2J) that generate the specific 
3d Lissajous curves m, a clear connection with multivariate polynomial 
approximation comes out, that could be useful in the corresponding data 
processing and analysis. 


Remark 3 (Clenshaw- Curtis type cubature). The availability of an hyper¬ 
interpolation operator with respect to a given density function (here the 
trivariate Chebyshev density) allows us to easily construct algebraic cuba¬ 
ture formulas for other densities, generalizing the Clenshaw-Curtis quadra¬ 
ture approach (cf., e.g., [26]). Indeed, if the “moments” 


= / 4>i,j,k( x ) £,(x)dx , i,j,k> 0 , i + j + k<n (20) 

are known, where £ E L\{{— 1, l) 3 ), as shown in [31] we can construct by 
(1131) the cubature formula 


f ^U n f{x)£(x)dx = Y C i, 3 , 
•'I -1 ’ 1 ] 3 0<i+j+k<n 


f~ 777 


= Y W sf{£n(0 S )) , W S =W S Y , (21) 

s=0 0 <i+j-\-k<n 

which is exact for all polynomials in P^. The resulting weights {Ws} are not 
all positive, in general, but if £/w E L 2 ((—1, l) 3 ), which is true for example 
for the Lebesgue measure £(x) = 1, it can be proved that 


lim Y \W 8 \ = / 


£ 0*0 

w(x) 


dx 


( 22 ) 


thus ensuring convergence and stability of the cubature formula; cf. m- 
We stress that these Clenshaw-Curtis type cubature formulas are based 
on Lissajous sampling (see Remark 2), and by Proposition 2 can be con¬ 
structed by a single 1-dimensional discrete Chebyshev transform along the 
Lissajous curve (i.e., by a single 1-dimensional FFT). 
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Remark 4 (Weakly Admissible Meshes and Discrete Extremal Sets). In 
the recent literature on multivariate polynomial approximation, the notion 
of “Weakly Admissible Mesh” has emerged as a basic tool, from both the 
theoretical and the computational point of view; cf., e.g., mmm and the 
references therein. 

We recall that a Weakly Admissible Mesh (WAM) is a sequence of finite 
subsets of a multidimensional (polynomial-determining) compact set, say 
A n C K C M. d (or C d ), which are norming sets for total-degree polynomial 
subspaces, 

\\p\\oc,K < C(An) IIpIIoo,^ > Vp G p£ , (23) 

where both C(A n ) and card(M„) increase at most polynomially with n. 
Here, denotes the space of d-variate polynomials of degree not exceeding 
n, and ||/||oo,x the sup-norm of a function / bounded on the (discrete or 
continuous) set X. Observe that necessarily card(M n ) > dim(P^). 

Among their properties, we quote that WAMs are preserved by affine 
transformations, can be constructed incrementally by finite union and prod¬ 
uct, and are “stable” under small perturbations [29]. It has been shown in 
the seminal paper HQ that WAMs are nearly optimal for polynomial least- 
squares approximation in the uniform norm. Moreover, the interpolation 
Lebesgue constant of Fekete-like extremal sets extracted from such meshes, 
say T n (that are points maximizing the Vandermonde determinant on A n ) , 
has the bound 

A(T n ) < dim(P^) C(A n ) . (24) 

Now, the Chebyshev lattices 

An = {£ n (e s ), s = 0,...,/x} (25) 

in (ED-®, form a WAM for K = [—1, l] 3 , with C(A n ) = 0((logn) 3 ). In 
fact, the corresponding hyperinterpolation operator 7i n being a projection 
on IP 3 , we get by (fTUD 

lblloo,[-l,l]3 = PnP||oo,[-l,l]3 < ll^nll ||p||oo,A. = C , ((logn) 3 ) \\p\\oo,An ■ 

( 26) 

In the next Section, we shall use the fact that Fekete-like extremal sets 
extracted from A n = {£ n (8s), s = 0,.. .,//} provide a Lissajous sampling 
approach to trivariate polynomial interpolation. 


4 Implementation and numerical examples 

4.1 Hyperinterpolation by Lissajous sampling 

In view of Proposition 2, hyperinterpolation on the Lissajous curve can be 
implemented by a single 1-dimensional Discrete Chebyshev Transform, i.e., 
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by a single 1-dimensional FFT. We shall concentrate on sampling at the 
Chebyshev-Lobatto points, since in this case we can conveniently resort to 
the powerful Chebfun package (cf. |21j). Sampling at the Chebyshev zeros 
can be treated in a similar way. 

Indeed, in view of a well-known discrete orthogonality property of the 
Chebyshev polynomials, the interpolation polynomial of a function g at the 
Chebyshev-Lobatto points can be written as 

n 

7i>(*) = ^ c m T m (t) (27) 

m =0 


where 

2 M 

c m = - " T m (r s ) g{r s ) , m = 1,..., p - 1 , 

/ J .s=0 

1 ^ 

c m = - T m (r s ) g(r s ) , m = 0,/r, (28) 

M s=0 

the double prime indicating that the Hrst and the last terms of the sum have 
to be halved (cf., e.g., [26, §6.3.2]). 

Applying this interpolation formula to g(t) = f(T an (t),Tb n (t), T Cn {t )) 
and comparing with the discrete Chebyshev expansion coefficients (1191) . we 
obtain by easy calculations 


7m 

'T'm 


f Cm , m = 1, 1 

7TC m , m = 0,+ 


(29) 


i.e., the 3-dimensional hyperinterpolation coefficients (11811 can be computed 
by the {c m } and PUI) . 

The coefficients of Chebyshev-Lobatto interpolation (1281) are at the core 
of the Chebfun package, cf. [L, 133]. A single call to the Chebfun basic 
function chebfun on f(T an (t),Tb n (t),T Cn (t)), truncated at the (/r + l)th- 
term, produces all the relevant coefficients {c m } in an extremely fast and 
stable way. 

For example, by the Matlab code m we can compute in about 1 second 
the p = |n 3 + |n 2 + n + 2 = 765102 coefficients for n = 100 with functions 
such as 


/i(s) = exp(—c||a;|||) , c > 0 , f 2 (x) = \\x\\% , /3 > 0 , (30) 


from which we get by (1131) the (n + l)(n + 2)(n + 3)/6 = 176851 coefficients 
of trivariate hyperinterpolation at degree n = 100. All the numerical tests 
have been made by Chebfun 5.1, in Matlab 7.7.0 with an Athlon 64 X2 Dual 
Core 4400+ 2.40GHz processor. 
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For the purpose of illustration, in Figure 1 we show the relative errors 
(in the Euclidean norm on a suitable control grid) for two polynomials of 
degree 10 and 20, respectively, and for the test functions f\ and fi in (1301) . 
Observe the Gaussian f\ is analytic, with variation rate determined by the 
parameter c, whereas the power function fi has finite regularity, determined 
by the parameter /?. 

Notice that the error decreases with the degree to a certain threshold 
above machine precision and thereafter does not improve. This is likely due 
to the fact that we require the summation of a large number of terms, for 
which a non-negligible error is to be expected. For practical applications 
this is of little import. 

In Figures 2 and 3 one can see the Chebyshev lattice on the Lissajous 
curve for polynomial degree n = 5. 


4.2 Interpolation by Lissajous sampling 

Concerning polynomial interpolation in the cube by sampling on the Lis¬ 
sajous curve, we resort to the approximate versions of Fekete points (points 
that maximize the absolute value of the Vandermonde determinant) studied 
in several recent papers U 0 [32]. By (1241) . it makes sense to start from 
a WAM, namely the Chebyshev lattice A n in (12511 . by the corresponding 
Vandermonde-like matrix 

V = V{A n -, 0) € R MxN , M = card(A0 = /t + 1 , N = dim(P^) , (31) 

(cf. (ED-© for the definition of p), where 

4> = Wi,j,k} , = T i (x 1 )T j (x 2 )T k (x 3 ) , 0 <i + j + k<n, 

is the total-degree trivariate Chebyshev orthogonal basis, suitably ordered 
(we adopt the graded lexicographical ordering , that is the lexicographical 
ordering within each subset of triples (i, j. k ) such that i + j + k = r,r = 
0,...,n). The (p, q) entry of V is the q-th element of the ordered basis 
computed in the p-th element of the nodal array. We recall that the choice 
of the Chebyshev orthogonal basis allows to avoid the extreme ill-conditiong 
of Vandermonde matrices in the standard monomial basis. 

The problem of selecting a N x N square submatrix with maximal deter¬ 
minant from a given M x N rectangular matrix is known to be NP-hard ESI, 
but can be solved in an approximate way by two simple greedy algorithms, 
that are fully described and analyzed in [6|. These algorithms produce two 
interpolation nodal sets, called discrete extremal sets. 

The first, that computes the so-called Approximate Fekete Points (AFP), 
tries to maximize iteratively submatrix volumes until a maximal volume 
N x N submatrix of V is obtained, and can be based on the famous QR 
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Figure 1: Top: Hyperinterpolation errors for the trivariate polynomials 
11 x 11 2 k with k = 5 (diamonds) and k = 10 (triangles), and for the trivariate 
function fi with c = 1 (squares) and c = 5 (circles). Bottom: Hyperinter¬ 
polation errors for the trivariate function /2 with /3 = 5 (squares) and f3 = 3 
(circles). 


factorization with column pivoting [8], applied to V f (that in Matlab is im¬ 
plemented by the matrix left division or backslash operator, cf. f27]). See 
m for the notion of volume generated by a set of vectors, which gener¬ 
alizes the geometric concept related to parallelograms and parallelepipeds 
(the volume and determinant notions coincide on a square matrix). 

The second, that computes the so-called Discrete Leja Points (DLP), 
tries to maximize iteratively submatrix determinants, and is based simply 
on Gaussian elimination with row pivoting applied to the Vandermonde-like 
matrix V. 

Denoting by A the M x 2 array of the WAM nodal coordinates, the 
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corresponding computational steps, written in a Matlab-like style, are 

w = V\v; s = f ind(u; / 0); Jyf FP = A(s,:); (32) 

for AFP, where v is any nonzero JV-dimensional vector, and 

[L, U, a] = LU(V, “vector”); s = <x( 1 : N)\ T PLP = A(s,:); (33) 

for DLP. In (13311 , we refer to the Matlab version of the LU factorization 
that produces a row permutation vector. In both algorithms, we eventually 
select an index subset s = (si,..., sjv), that extracts a Fekete-like discrete 
extremal set J- n of the cube from the WAM A n . 

Once the underlying extraction WAM has been fixed, differently from the 
continuum Fekete points, Approximate Fekete Points depend on the choice 
of the basis, and Discrete Leja Points depend also on its order. An important 
feature is that Discrete Leja Points form a sequence , i.e., if the polynomial 
basis is such that its first N r = dim(Pj?) elements span Pf, 1 < r < n (as it 
happens with the graded lexicographical ordering of the Chebyshev basis), 
then the first N r Discrete Leja Points are a unisolvent set for interpolation 
in Ff. 

Under the latter assumption for Discrete Leja Points, the two families 
of discrete extremal sets share the same asymptotic behavior, which by a 
recent deep result in pluripotential theory, cf. [2], is exactly that of the 
continuum Fekete points: the corresponding uniform discrete probability 
measures converge weakly to the pluripotential theoretic equilibrium measure 
of the underlying compact set, cf. .4, 6j. In the present case of the cube, 
such a measure is the product Chebyshev measure with scaled density 
w(x)/ir 3 . 

We give now some numerical examples, that can be reproduced by the 
Matlab package [18] . First, in Figures 2-3 we show the Approximate Fekete 
Points extracted from the Chebyshev lattice on the Lissajous curve for degree 
n = 5. In Figure 4 we display the numerically evaluated Lebesgue constants 
of the Approximate Fekete Points and Discrete Leja Points for degree n = 
1,2,... ,30. For both the nodal families, the Lebesgue constant turns out 
to be much lower than the upper bound (1241) . and even lower than N = 
dirn(P^), a theoretical upper bound for the continuum Fekete points. In 
particular, the Lebesgue constant of Approximate Fekete Points seems to 
increase quadratically with respect to the degree, at least in the given degree 
range. 

Finally, In Figure 5 we show the relative interpolation errors for the two 
test functions f\ and /2 of Figure 1. Since the Discrete Leja Points form a 
sequence, as discussed above, we have computed them once and for all for 
degree n = 30, and then used the nested subsequences with N r = dirn(P(?) 
elements for interpolation at degree r = 1,... ,30. The corresponding file 
of nodal coordinates can be downloaded from [18] . The relevant indexes 
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(si, S 2 , ■ ■ ■, sn 30 ) corresponding to the extraction of the Discrete Leja Points 
from the Chebyshev lattice d25l) -(l9l) at degree 30, could be used in appli¬ 
cations, such as MPI |25], where a trivariate function is not known or 
computable everywhere, but can be sampled just by travelling along the 
Lissajous curve. 



Figure 2: The Chebyshev lattice (circles) and the extracted Approximate 
Fekete Points (asterisks), on the Lissajous curve for polynomial degree n = 5. 
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5 Appendix 

Proof of Theorem 1. We prove the theorem for n even, the proof being 
similar in the odd case. Let be m = n/2, n even, so that 

( a n , b n , c n ) = (3m 2 + m, 3m 2 + 2m, 3m 2 + 3m + 1) . 

We assume that 

i 2 + j 2 + k 2 >0. 

First case. We show that it is not possible to have 


ia = jb + kc 


for i + j + k < 4m (= 2 n). Now, ia = jb + kc becomes i(3m 2 + m) + 
j (3m 2 + 2m) + k{3m 2 + 3m) + k. Since m divides 3m 2 + m, 3m 2 + 2m and 
3?n 2 + 3?n, we must have that m divides k, i.e., k = am, a > 0. Since 
k < 4m, 0 < a < 4. 

Hence we must have 

i(3m 2 + m) = j(3m 2 + 2m) + am(3m 2 + 3m + 1) 
that is, dividing by m, 

i(3m + 1) = j(3m + 2) + cc(3m 2 + 3m + 1) , 
which is equivalent to 

i((3m + 2) — 1) = j(3m + 2) + a((3m + 2)m + (m + 1)) 


and to 


(3m + 2)(i — j — ma) = i + a(rn + 1) . 


The latter implies that 


i + a(m + 1) = /3(3m + 2) 
for some integer f3 > 0, i.e., 


i = /3(3m + 2) — a(m + 1) 


(actually /3 = i — j — ma). 
From 


f3 = i — j — ma 
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we have 


j = i — ma — (3 = /3 (3m + 2) — a(m + 1) — ma — (3 

i.e., 

j = (3(3m + 1) — a(2m + 1) 

(which must be > 0). It follows that 

z + j + k = (3(3m + 2) — a(m + 1) + (3(3m + 1) — a(2 m + 1) + am , 

i.e., 

i + j + k = /3(6m + 3) — a{2m + 2) . 

We now consider two possibilities for a: 

1) a = 0. In this case 

i = (3(3m + 2) , j = [3(3m + 1) , k = 0 

and i + j + k = [3(6m + 3). Now, /9 / 0 otherwise i = j = k = 0. 
Hence 

z + j + k > 1 (6m + 3) > 4m 
violating the constraint on i + j + /c. 

2) a > 1 (and a < 4). In this case (3 > 1, for otherwise i, j < 0. More 
precisely, since 

j = j3(3m + 1) — a(2m + 1) = (3/3 — 2a)m — a > 0 
we must have 3/3 — 2a > 1. Hence 

z + j + fc = /3(6m + 3) — a(2m + 2) = m(6/3 — 2a) + 3/3 — 2a 

= m(3/3 — 2a + 3/3) + 3/3 — 2a > m(l + 3) + 1 = 4m + 1 > 4m 
which again violates the constraint on i + j + k. 

Second case. It is not possible that 

jb = ia + kc 

for i + j + k < 4 m (= 2n). In this case, ia = + kc becomes i(3m 2 + m) = 

j(3m 2 + 2m) + k(3m 2 + 3m) + fc. Since m divides 3m 2 + m, 3m 2 + 2?n and 
3 ?b 2 + 3m, we must have that m divides k, i.e., k = am, a > 0. Since 
k > 4m, 0 < a < 4. 

Hence we must have 

j(3m 2 + 2m) = i(3m 2 + m) + am(3m 2 + 3m + 1) 
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and dividing by m 

j(3m + 2) = i(3m + 1) + a(3m 2 + 3 m + 1) 
which implies that 

j(3m + 1) + j = i(3m + 1) + a(m(3m + 1) + 2m + 1) 

and also 

j — a(2m + 1) = (i — j + am)(3m + 1) . 

Let /3 = i — j + am (which a priori could be < 0) so that 

j — a((2m + 1) = (3(3m + 1) 

which is equivalent to 


j = j3{3m + 1) + a{2m + 1) , 


and 

i = f3 + j — am = f3 + (/ 3{3m + 1) + a(2m + 1)) — am , 

i.e., 

i = /3(3 m + 2) + a(m + 1) . 

Hence 


% + j + k = f5{3m + 2) + a(m + 1) + j3(3m + 1) 
+ a{2m + 1) + am 
= (3 (6m + 3) + a (4 m + 2) 

= m( 6/3 + 4a) + 3/3 + 2a 
= (3/3 + 2a) (2m + 1) . 

For 0 < i + j + k < 4m, the only possibility is 

3/3 + 2a = 1 . 

For 0 < a < 4, the only integer solution for f3 is 

a = 2 , /? = -!. 


However, in this case, 

i = (3(3m + 2) + a(m + 1) = —(3m + 2) + 2 (m + 1) = — m < 0 
which is not allowed. 

Third case. It is not possible that 

kc = ia + jb 
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for i+j+k < 4 m (= 2n). In this case, kc = ia+jb becomes k(3m 2 +3m)+k = 
i(3m 2 +m)+j(3m 2 +2m). Since m divides 3m 2 +m, 3m 2 +2m and 3m 2 +3m, 
we must have again that m divides k, i.e., k = am, a > 0. Since k > 4m, 
0 < a < 4. 

Hence 


am(3m 2 + 3m + 1) = i{3m 2 + m) + j(3m 2 + 2m) . 
Dividing by m we obtain 

a(3m 2 + 3m + 1) = i(3m + 1) + j(3m + 2) 

or equivalently 

a(m(3m + 2) + m + 1) = i(3m + 2 — 1) + j(3m + 2) 


and 


Let )3 


i + a(m + 1) = (3m + 2 )(—am + ?' + j) . 
—am + i + j. Then 

i + a(m + 1) = /3(3m + 2) 


which implies that 

i = (3(3m + 2) — a(m + 1) = m(3/3 — a) + ( 2/3 — a) . 

Note that i > 0 implies j3 > 0 (since a > 0). Further 
j = /3 + am — i = (3 +am — (/3(3m + 2) — a(m +1)) = a(2m +1) — /3(3?n +1) , 


i.e., 

and 


j = m{2a — 3/3) + (a — (3) 


i + j + k = f3(3m + 2) — a(m +1) + a(2m +1) —/3(3m+ 1) + am = /3-\-2am . 
If a = 0, then 

i = /3(3m + 2) , j = —f3(3m + 1) , k = 0 

which is not allowed as j > (1 (and (3 > 0). 

If a = 3,4 

i + j+ k = /3 + 2 am > 6 m > 4m 

which also contradicts the constraints on i+j + k. 

If a = 2, 

i + j + k = (3 + 4m > 4 m 
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unless P = 0. However, in this case 


i = —2 (m + 1) < 0 


and so a = 2 is not possible. 

The only remaining possibility is a = 1. In this case 

i = (3(3m + 2) — (m + 1) , j = (2 m + 1) — /3(3 m + 1) , k = m . 

But j > 0 is equivalent to 2 m + 1 > (3(3m + 1), i.e., 

2m + 1 . 1 

P < - < 1 , for m > 1 

3 m + 1 

and so P = 0 (as P is an integer). But then 

i = — (m + 1) < 0 


which is not possible. 

Counterexample. Let 


i = 2m + 1 , j = m , k = m . 

Then i + j + k = 4 m + 1 and it is elementary to check that ia — jb — kc = 
0. Hence, 4m = 2 n is the maximal value for which the property in the 
statement of Theorem 1 is satisfied. □ 
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